Draft version March 3, 2009 

Preprint typeset using 1^1^^ style cniulatcapj v. 10/09/06 



TIME-DEPENDENT MODELING OF RADIATIVE PROCESSES IN HOT MAGNETIZED PLASMAS 

Indrek Vurm^ and Juri Poutanen 

Astronomy Division, Department of Physical Sciences, P.O.Box 3000, 90014 University of Oulu, Finland; indrek.vurm@oulu.fi, 

juri.poutanen@oulu.fi 
Draft version March 3, 2009 

ABSTRACT 

Numerical simulations of radiative processes in magnetized compact sources such as hot accretion 
disks around black holes, relativistic jets in active galaxies and gamma-ray bursts are complicated 
^^ ' because the particle and photon distributions span many orders of magnitude in energy, they also 

strongly depend on each other, the radiative processes behave significantly differently depending on the 
energy regime, and finally due to the enormous difference in the time-scales of the processes. We have 
developed a novel computer code for the time-dependent simulations that overcomes these problems. 
The processes taken into account are Compton scattering, electron-positron pair production and 
annihilation, Coulomb scattering as well as synchrotron emission and absorption. No approximation 
has been made on the corresponding rates. For the first time, we solve coupled integro-differential 
kinetic equations for photons and electrons/positrons without any limitations on the photon and 
lepton energies. A numerical scheme is proposed to guarantee energy conservation when dealing 
with synchrotron processes in electron and photon equations. We apply the code to model non- 
f^ , thermal pair cascades in the blackbody radiation field, to study the synchrotron self-absorption as 

Oh' particle thermalization mechanism, and to simulate time evolution of stochastically heated pairs and 

corresponding synchrotron self-Compton photon spectra which might be responsible for the prompt 
^ , emission of gamma-ray bursts. Good agreement with previous works is found in the parameter regimes 

^ ' where comparison is feasible, with the differences attributable to our improved treatment of the 

j^ I microphysics. 

Subject headings: accretion, accretion disks — galaxies: active — gamma rays: bursts — gamma rays: 
_- 1 theory — radiation mechanisms: non-thermal — X-rays: binaries 
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t:^ ; 1. INTRODUCTION 

^^ ' Spectral energy distributions of a number of compact, magnetized, high-energy sources such as relativistic jets from 
^^ , active galaxies, gamma-ray bursts, black hole accretion disk-coronae are stron gly affected and shaped by Compton 
t — ■ scattering, synchrotron radiation and electron-positron pair production (see e.g. lGierlihski et aLlllQQgtrGhiselhni et al.l 
^^ ' I2OO2I : IZdziarski fc Gierlihskill2004 IStern fc PoutanenI 12004 ). Understanding the physical conditions in these sources 
00 , requires detailed modeling of the interactions between the particles and photons, which is not an easy task. The basic 
^^ ' problem is that we cannot compute radiative processes from a given a priori lepton distribution (e.g. Maxwellian 
or a power-law), because it depends strongly on the radiation field, which in its turn is determined by the particle 
distribution. Another problem is that the time-scales for various processes differ by orders of magnitude. The energy 
range of particles and photons responsible for the emission also spans many orders of magnitude, with different 
|h processes making dominant contributions to the emergent spectrum in different bands. One of the main difficulties 
in calculating radiative processes over a wide range of energies is that a particular radiative process may behave 
significantly differently depending on the energy regime, the most well-known example of such processes being also 
the most important in relativistic plasma, namely Compton scattering. Depending on the energies of the interacting 
particles, an electron or a photon can lose (or gain) a significant or negligible fraction of its initial energy in one 
scattering. The former case has to be accounted for by the integral scattering terms in the kinetic equations, while 
the latter necessitates the Fokker-Planck treatment. 

The treatment of radiative processes in re lativistic pla smas has been the subject of several works. There are two 
basic approach es: Monte Carlo methods (e . g.lStern et aLlll9 95: Pilla & Shaham 1997) and solving the relevant kinetic 
equations (e.g. lLightman fcZdziarskilll987l : iCoppilll99a 11999, : .Navakshin fc MeUai.l998l : lPe'er fc Waxmanll2005D . Both 
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have their own advantages and disadvantages. Monte Carlo treatment makes it easy to take into account radiative 
transfer effects, on the other hand it usually suffers from poor photon statistics at high energies. Another problem 
can arise at very low energies, where the optical thickness to synchrotron absorption can be enormous. In the kinetic 
theory approach photon statistics is not an issue, the sole difficulty lies in solving the relevant integro-differential 
equations. In this work we have chosen to follow the second approach. 

Due to the difficulties in solving the exact Boltzmann equations of the kinetic theory, different simplifying ap- 
proximations have been made in earlier works. They fall in three basic categories: ad hoc assumptions about the 
particle energy distributions, approximate treatment of different physical processes, and simplified treatment of ra- 
diative transport. Various approximations invoked to simplify the treatment of radiative processes at the same time 
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limit the range of their apphcabihty. One commonl y employed approxima tion concerns Compton scattering, which is 
assumed to take place in the Thomson regime (e.g. iGhisellini et al.l[l998L hereafter GHS98) and is accounted for by 
a simple cooling term in the electron equation. This sets two restrictions that the photon energy in the electron rest 
frame is smaller than the electron rest mass and the average photon energy is much lower than the electron kinetic 
energy. Otherwise all photons would not contribute to electron cooling, the higher energy ones being downscattered via 
Compton scattering. This means that one is unable to treat cases when Comptonization approaches saturation, which 
may be relevant at high compactnesses, and to study electron heating by external radiation. Other works account 
for Klein-Nishina corrections to the electron cooling rate, but still neglect the diffusive nature of the process when 
electron and photon energies are comparable (Coppi 1992; Moderski et al.ll2005t iPe'er fc Waxmanll2005l ). which works 
towards establishing an equilibrium Maxwellian distribution. Another useful approximation, when the integral terms 
desc ribing Compton scattering are accoun ted for, is to consider ultrarelativistic electrons and very low energy photons 
(e.g. lZdziarskil[l988ilModerski et al.ll2005( l. This, however, becomes increasingly inaccurate when the electrons cool to 
sufficiently low energies. 

The cyclo-synchrotron process also exhibits qualitatively different behavior depending on the energy of the radiat- 
ing particles. If the emitting particles are relativistic, the emission spectrum is smooth and can span several orders 
of magnitude in energy, while in the nonrelativistic case the energy is radiated at discrete cyclotron harmonics and 
most of this radiation might be strongly self-absorbed. In the first case, the radiating particle (electron or positron) 
mostly loses its energy in a continuous fashion, while in the second case it can gain energy by absorbing the cyclo- 
synchrotron photons emitted by other particles. T his process is a dominant particle thermalization mechanism in 
compact magnetized sources (IGhisellini et al.lll988f ). Its proper account requires accurate emissivities in the transrel- 
ativistic regime, because electron thermalization usually takes place at mildl y relativistic energ ies. Some codes for 
computing radiative processes in relativistic plasma (e.g. eqpair described in ICoppilll992l [l999l ) neglect this process 



completely as the electrons are as sumed to be thermal at low energies or account for thermalization by Coulomb 
collisions only (jNavakshin &: Melial [l998l. In other approaches synchrotron thermalization is computed (GHS98), but 
Compton scattering is then considered only approximately. 

Owing to the fact that proper treatment of transport processes for all types of particles would make the task 
prohibitively difficult, and partly due to our ignoran ce of the exact geometry of the probleni , it is rather common 
practice to neglect radiative transport altogether (e.g. iLightman fc Zdziarskilll987l : ICoppi|[l99^ and assume spatially 
homogeneous and isotropic particle distributions. In this case particle and photon loss from the system is modeled in 
terms of simple escape probabilities. We also follow this approach here. 

In this paper we introduce and describe in details a numerical code that can deal with Compton scattering, syn- 
chrotron emission and absorption, electron-positron pair production and annihilation without limitations on the en- 
ergies of the photons and electrons/positrons. We solve coupled integro-differential kinetic equations describing time 
evolution of the photons and lepton distributions simultaneously. When necessary, the Fokker-Planck differential 
terms are substituted instead of the integral terms with coefficients computed exactly from the moments of the inte- 
gral equation. Particle thermalization by synchrotron self-absorption, Coulomb (M0ller) scattering as well as Compton 
scattering is considered. Extreme caution is taken when dealing with synchrotron self-absorption, because of cancel- 
lation of large, almost equal terms, which can result in inaccuracies and huge energy sinks. Numerical simulations 
show that our code conserves energy with about 1% accuracy. We present an extensive testing of the code using some 
problems described previously in the literature. Processes involving bremsstrahlung can be easily added to the code, 
while for the conditions considered in the paper, they are not important. 

2. THE KINETIC EQUATIONS 

We are considering a region of relativistic plasma of charged particles (electrons and positrons, which we call "elec- 
trons" below if the relevant processes, e.g. Compton scattering and synchrotron, operate identically on both types of 
particles) permeated by radiation and tangled magnetic fields. We study the evolution of lepton and photon distri- 
butions by solving the time-dependent coupled kinetic equations accounting for synchrotron emission and absorption, 
Compton scattering. Coulomb scattering, and electron-positron pair production and annihilation. We make a simpli- 
fying assumption of homogeneity and isotropy of the particle distributions. We assume that energy is transferred to 
electrons by some unspecified mechanism, which is manifested as either injection of high-energy electrons or diffusive 
acceleration within the active region. The escape of radiation (and also electrons) from the region is modeled by a 
simple escape probability formalism. 

2.1. Distribution functions 

The dimensionless four-momentum of a photon is x = {x,x} = a;{l,ti;}, where u) is the unit vector in the photon 
propagation direction and x = hv/rac(? . The photon distribution can be described by the occupation number nph or 
by the photon number density per linear and logarithmic interval of photon energy: 

A^ph = / Nph(x) dx = Uphix) dlnx — —3- / d'^uj / np]^{x) x^ dx, ( 2-1) 

where Ac = h/mcC is the Compton wavelength. Functions A^ph(x) and fiph are used in general forms of kinetic 
equations and r7,ph(x) is convenient for numerical work. 
The dimensionless electron (positron) four-momentum is defined as p = {7,^} = {7,^11} = 7(1, /3J1}, where H is 

the unit vector in the electron propagation direction, 7, /3, and p — (3j = \/j'^ ~ 1 are the electron Lorentz factor. 
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dimensionless velocity and momentum, respectively. We can use subscripts + and — to distinguish between positrons 
and electrons. The electron/positron distributions can be defined in a number of alternative ways (normalized to their 
number density): 

N± ^ J N±{j)dj ^ J n±{p) dlnp ^ ^ J d'n J n±{p) p^dp. (2-2) 

The occupation number n±(p) and the density per unit Lorentz factor are useful quantities used in general kinetic 
equations, while the electron density per logarithmic momentum interval n±(p) is more appropriate for numerical 
work. For the processes, where the distinction between electrons and positrons is unnecessary, we use the sum of the 
distributions, for example, ric = n_ + n+. 

2.2. General form of the kinetic equations 

The relativistic kinetic equation (RKE) describing the evolution of the occupation number ni(pi) of sp ecies 1 
(electron or photon) as a result of binary collisions can be written in the covariant form (|de Groot et aLNlOSClD 

f d^p2 d^p^ d^p4 
p • Vni(pi) = / 6{p + p^ - p^ - p,)Wi2^34:[n3{P3)ni{p4^) ~ ni{p^)n2{p2)] , ( 2-3) 

where V = {d/cdt, V} is the four-gradient, e^ is the zeroth component of the corresponding four-momentum, and 
M^i2^34 = W^34^i2 is a Lorentz scalar transition rate, which possesses the obvious symmetry. In this equation, the 
non-linear terms related to fermion degeneracy and induced photon scattering are omitted. As it stands, the right-hand 
side of equation (|2-3p accounts for the rate of one particular process. To determine the full evolution of fii we should 
therefore sum up the coUisional integrals accounting for all relevant processes. 

In the frame where the particle distributions are isotropic (we call this frame E), the kinetic equation can be 
represented in the form (skipping subscript 1): 



dfilv) 1 (9 r \ 1 DniT)] 

' ^^~ i eeph(p) - -— [D{e)epn{p)] - 



dt p^dp{ "^ ''^' 2de' ' ' "^ '^'^j Dt 



( 2-4) 
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where the momentum derivative term accounts for continuous energy gain/loss processes, while the right hand side 
contains all discontinuous processes such as scattering, emission, absorption and escape. The quantities e and -D(e) 
account for systematic particle heating/cooling and diffusion in energy space, respectively. Both are generally energy- 
dependent for the processes we are considering here. For the following discussion it is convenient to decompose the 
kinetic equations in terms of the contributions from different physical processes as 



dt 

dn±{p) 

dt 



n-ph 


syn(a;) + 


"ph 


cs(a;) 


+ »^pii 


pp(a^) 


'-CSC 


Qph, 


"i, 


3yn(p) + 


n±,c 


,{P) + 


'^i.PF 


b) + 


n±,Coni{p) - 


n±{p) 



2-5) 



where syn, cs, pp and Coul stand for synchrotron, Compton scattering, pair production (and annihilation), and 
Coulomb scattering, respectively. The terms describing physical processes can contain both differential and integral 
parts, depending on the nature of the process and the way we find most convenient to treat it. Thus the equation for 
photons has the form: 



a?iph(a;) _ d 



dt d\nx 



^ph(a;)r7,ph(a;) - Sph(x)^-2 

omx 



ifph(a;,xi)7iph(a;i) dlna:i V S-pi,. ( 2-7) 

tph 



Here the differential term is responsible for Compton scattering in diffusion approximation, while the integral term with 
kernel -ftTph describes scattering that can be resolved on the grid. The sink term ex 1/iph describes photon absorption 
(by synchrotron and pair-production) and scattering as well as the escape, while S'ph gives the contribution from pair 
annihilation, synchrotron emission and other (e.g. blackbody) photon injections. 
Similarly for electrons and positrons, we write 



dn± [p) d 



dt 9 In p 



Mp)n^{p)-B,{pf''^^^'^ 



d\np 



/ Kc{p,pi)n±{pi) dlnpi 



n±{p) 



+ S±, 



where coefficients A^ and B^ describe electron cooling, heating and diffusion as a result of synchrotron emission and 
absorption, Compton scattering in Thomson limit, Coulomb scattering as well as possible diffusive particle acceleration. 
The integral term with kernel K^ describes Compton scattering in Klein-Nishina limit into the bin and the sink term 
ex l/i± gives the scattering from the bin as well as the electron escape and pair annihilation. The source term S± 
contains pair production as well as a possible electron injection term. 

2.3. Escape probability formalism 

As we are studying radiative processes in a simple one-zone framework neglecting the radiative transport effects, we 
must include an escape term in equation (|2-5p to allow for the fact that photons can leave the emission region of finite 
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size R and produce the radiation flux that is actually observed. The typical escape timescale is usually estimated from 
random walk arguments resulting in tcsc ^ R{1 + Tsc)/c, where Tsc is the scattering opacity. Such form accounts for the 
fact that if multiple scatterings are important (tsc > 1), photons have to 'diffuse' out of the medium and the escape 
time is prolonged by a factor Tsc- However, it does not account for the fact that if the medium is absorptive, a typical 
photon cannot diffuse further than the thcrmalization length l^, — [aa(cKa + Qfsc)]"^ before it is destroyed (oa and 
asc are extinction coefficients due to absorption and scattering, respectively). To incorporate both effects, we employ 
the solution of a simple radiative diffusion problem in a sphere of radius R with constant emissivity, absorptivity 
and monochromatic scattering. The escape timescale is estimated by comparing the emergent flux to the radiation 
density inside the source. While clearly an oversimplification, such estimation nevertheless has the desired properties 
mentioned above. 

Defining the effective optical thickness of the medium as r» — ^3ra(Ta + Tsc), where t^ = a^R and Tsc = ctscR are 
optical thicknesses due to absorption and scattering, respectively, we find 



_2R j V3 



T* ( 1 — C 



■) 



r*(l + e-2^-)-(l-e-2^*) 



3 



(2-9) 
reduces 



where A = asc/icta + asc) is the single-scattering albedo. If the medium is translucent (r, <^ 1), equation 
to a more familiar form 

op/ o 

tcsc = — I 1 + —Tsc ] ■ ( 2-10) 
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In our simulations, aa includes cyclo-synchrotron absorption and photon-photon pair production, and asc is the 
extinction coefficient for Compton scattering. 

2.4. Compton scattering 

2.4.1. Compton scattering of photons 

The explicitly covariant form o f RKE for Compton scattering of photons ignoring non-linear terms is (|Pomraninsl 
119731 : iNagirner fc Poutanenlfl993 . hereafter NP94) 

d(p +Xi - p - x) F [nph{xi)nc{Pi) - np\-,{x)nc{p)\ , 

^c J 7 71 ^1 ~^ 

where r^ is the classical electron radius, F is the Klcin-Nishina reaction rate ([Berestetskii et al]|1982D 

2 



r^ 2 /■ 
x-Yriph{x) = y^ / ■ 



( 2-11) 
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and S, = p ■ Xi = p ■ X and £,i = p. ■ x= p ■ x^ are the scalar products of four- vectors. 

We assume the existence of a reference frame where the particle and photon distributions are approximately homo- 
geneous and isotropic. Under the spacial homogeneity assumption we can write equation (|2-lip as 

Dn \i{x) 1 [ (Pxi 

E- — =-caTSo{x)NcnpY,{x)+caTNc- Rph{xi ^ x) npY,{xi). (2-13) 

^^ coll,cs ^ J ^1 

The scattering cross-section (in units of Thomson cross-section ctt) is given by 

3 2 1 



so{x) = 
and the redistribution function is 



167r \l.Nc X 



d^p (fpi (fxi ~ f . rpxf , ^ 

7io(p) F d[p^+x^-p-x) 



7 71 xi 



Rpb{xi -^x)= / ndPi) FS{p +XJ^~p- 

16^T X-l^Nc J 7 71 -1 _ 

For isotropic particle distributions in frame -E, equation (|2-13p can be written as 



Dhpi,{x) 



Dt 



47r 

'CaTloix) NcfipYiix) + cutNc — / xidxi Rp\,{x,xi) np\,{xi) 



( 2-14) 



( 2-15) 



( 2-16) 



where the redistribution function averaged over the cosine of the scattering angle ^ ~ x-xi/{xxi) = a;-u;i is expressed 
via an integral over the electron distribution (NP94): 

— 1 /"^ 3 2 f°^ — 

Rph{x,xi) = - Rph{xi -^ x) dfj, ^ — ^^ ^^ j Rph{x,xi,'^i)hc{pi) dji. (2-17) 



Here 



i?ph(a;,a;i,7i) = ^ Pi 



16 A^iVo Jj^i^x.xi) 



d u 

d?Q.i Suji F 5{p -\- XjY — p — x), 



( 2-18) 
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and the lower limit of the second integral in equation (|2-17p comes from the condition of energy and momentum 
conservation: 

([x~xi + {x + xi)y/TTT/xx^]/2 if |a;-xi| > 2a;a;i, r 9 10^ 

^*^'^'^''^\l + {x-xi + \x-xi\)/2 if |a;-a;i| <2x:ei. ^ ' 

The integrals in equation (|2-18p can be calculated analytically (jBrinkmannl ligS^I : NP94) to obtain a fully general 
expression fo£_Rph(a;, xi, 7i) valid in all regimes (see Appendix \B\i . This is an alternative form of the function derived 
bv [Jo^(ll96l . ^^ 

Since the total number of particles is conserved in Compton scattering, multiplying the rhs of equation (|2-16p by 
x"^ integrating over dx must give zero, implying a relation between the redistribution function and the extinction 
coefficient (NP94) 

so{x) = — / Rph{xi,x) xi dxi . ( 2-20) 



This can also be inferred directly from the definitions (|2-14p and (|2-15p . 

In the kinetic equation (|2-5p for the photon density nph{x) the Compton term is obtained by multiplying equation 
(12^161) by SnX^^x^. 

2.4.2. Compton scattering of electrons and positrons 

The description of Compton scattering for electrons and positrons is very similar to that for photons. In the linear 
approximation the RKE reads 

£• Vn±(p) = ^-3- / -S{p^ +£1 - p- x) F [fiph{xi)n±{pi) - fipb{x)n±{p)] . ( 2-21) 

z Aq j X X\ '7i 

Neglecting spatial gradients, equation (|2-2ip becomes 



Dn±{p) 



Dt 



coll,cs 

where the scattering cross-section for electrons is 



1 /■ d-^pi 

ccrTSo(p) A^phn-±(p) + ccrTiVph- / ^o(Pi ^p)".±(Pi), (2-22) 

1 J 7i 



3 2 1/" d^x d^xi d^pi _ , , ^ c-/ N / „ „„N 

and the redistribution function 

i?e(Pi-p) = ^^/^^nph(a.r)f%+a-p-£). (2-24) 

Making use of the isotropy of the problem, we can rewrite the kinetic equation in frame E for isotropic distribution 
no(p): 

Dn± (p) 



Dt 



Att f — 

= -caT so{p) Nph n±{p) + caTNph — / pidji Rc{p,pi) n±{pi), ( 2-25) 

coll,cs ' ^ 

where the electron redistribution function averaged over cosine of the electron scattering angle ^e is 

1 /"^ 3 2 f°° 

Rc{p,Pi)^7; Rc{Pi -^ p) dfie = 7^ .3 ., / i?o(7,7i,a;i)nph(xi) dxi, (2-26) 

where 

R,i-f,-fi,xi) = ^xi f -^d^ojid^n^FSip^+x.-p-x) (2-27) 

and 

a^47,7i) = [7-7i + b-Pi|]/2. (2-28) 

The relation between the redistribution function and the extinction coefficient is 

An f — 
so{p)^— J Re{Pi,p)Pi dii . (2-29) 

Not surprisingly, there turns out to be a relation between the quantities i?ph and Re (proved in Appendix[A|, namely 

ppiRe{'^,^i,xi) = xxiRpii{x,xi,^i), ( 2-30) 

together with the energy conservation condition a; + 7 = xi + 71. Through equations (|2-30p and (|2-18p we have a 
generally valid expression also for i?o(7,7i,a;i). 

In the kinetic equation (j2-6p for the electron and positron densities n±{jp) the Compton terms can be obtained by 
multiplying equation (|2-25p by ^nX^^p^ . 
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2.5. Photon-photon pair production and pair annihilation 

Th e electron RKE accounting for pair-production and annihilation processes can be written as ([Nagirner fc Loskutovl 
fl999l) 

r^ 2 /■ d^p+ d^Xi d^x 

p_ • Vn_(p_) = ^-3- / 6{p_+p -x^-x)F^^[fLph{xi)fLph{x)-il^{p_)fl+{p^)], (2-31) 

4 Aq j ^-j_ Xi X 

where we used subscripts =F to explicitly show the momenta and the occupation number of electrons and positrons. 
Assuming homogeneity, we get 

Dh-{p_) 



Dt 



coll,pp 

where the pair annihilation cross-section (in units of ctt) is given by 



= -c (IT Spa(p_) N+ n-{p_) + c fTT ^ph -f Jpp(P-), ( 2-32) 



and the pair production rate by 

2 



_ , , 3 2 If d^p+ d^xi d^x . , ^ „ f, ^ / „ „„N 

SpaP_) = ^T^Tg-TT / n+{p^)F^^6{p +p -xj^-x) (2-33 



3 f 2 Y I f d^p+ d^xi d^x ^ 



Jpp(P-) = ^TK-\ ^3 AT — / nph(a;)nph(a;i)F^^<5(p +P -x^^-x). (2-34) 

327r VAJiVph/ 7- J 7+ a;i a; _- _ + 

The relativistically invariant reaction rate Fj^ is (jBerestetskii et al1ll982t ) 



e . 6 , o ^1 1\ /I 1 



2 



^„ = ^ + J + 2^J + -J-^- + -J. (2-35) 

where ^ = p ' 2L — P ■ ^1 find ^1 = p ■ 2I1 = P ' HH- 
Assuming again isotropic particle distributions in frame E, we can write equation (|2-34p as 

jpp(p-) = Stt -3—— / np^{x)dx\ fiph{xi)dxi R^^{j_,x,xi), (2-36) 

VApiVph/ j-p- J^(L) J^iL) 

where we have defined 

i?^-y(7-,a:,a:i) ^ - xxip- / d^uj d^uji F^^ S{p _ + p - x^ - x). (2-37) 

2 (47rj^ J 7+ - - + 

The cross-section becomes 

2 Z""" 2 
Spa(p-) = 47r / p+(ip+ crpa(7+,7-) n+{p+), ( 2-38) 



'P---^--"A^iV, 



where 



o'pa(7-H,7-) = oTTZTI / d^f7+F^^<5(£_+£ -xj-x). (2-39) 



d^xi d^x 

(47]-)2 ^_'y^ J Xi X 

The treatment of positrons is identical if we switch the subscripts — and + in equations (l2-3ip - p-39p . 
The photon kinetic equation accounting for pair production/annihilation processes is 

^- V7iph(a;) = ^-3- / i 5{j) +p -x^ -^)i^^,^ [fi_(p_)fi+(p+) -nph(a;i)nph(a;)] . (2-40) 

2 AJ J xi 7- 1+ -- - + 

Neglecting the spatial derivatives in the left hand side, this becomes 



Dnp\^(x) 



A 



Dt 
where the pair-production cross-section is 



= -c CTT spp(a;) TVph nph(a;) +c<7tN-N+-^ jpa(a;), ( 2-41) 

coll.pp 



3 2 I f d^xi d^p- d^p+ ^ 



spp(a;) = —- -3— / fipY,{xi) F^^ 5{p +p -x^~x) (2-42) 

le-TT AJTVph X J xi 7- 7+ -- - + 



and the emissivity due to pair annihilation 



3 / 2 V 1 If d^xi d^p- d^p 



-1+ ~ 



Jpa(a;) = TTT TT TTITt / n-iP-) "+(P+) ^77 ^(P +P,-^i-x). ( 2-43) 

loTT \XqJ N-N+xJ Xi 7- 7+ -^ - + 

Notice that unlike the electron equation, the photon equation is nonlinear owing to the fact that the cross-section 
(|2-42l) depends explicitly on the photon distribution. 
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Under the isotropy assumption equations (|2-42p and (|2-43p in frame E become 



2 \^ 1 1 



■VC/ ---+- J7i^' ^7i^' 



jpaW = 67r ( ^ 1 jy j^ ^ I „,"+(P+)c^7+ I ,^^n_iP-)d-l-Ry^h^,x,Xi), (2-44) 



where we have to substitute cci = 7_ + 7+ — x from the energy conservation condition, and 
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,(a;) = An I x\dxi ap-p{x,xi) nph(a;i), ( 2-45) 

^c^Vph Ji/x 



app{x,xi) = -jjz^zzr I ~:. T"*" ^^^1 ^fi ^^P- +P+ ^^1 ~2)- ( 2-46) 



where 

4 (47r)2 xxi 7 7_ 7-f 

Explicit expressions for the rate R^^{'^-,x,xi) (derived bv ISvenssonlll982L see also iBoettcher fc Schlickeiserlll997l and 
iNagirner fc Loskutov"1999') and the cross-sections ov,a (7+, 7-), crpp(x,xi) as well as the lower integration limits in 
equations (|2-36p and (|2-44p are given in Appendix |D] 
The pair-production terms in equations (|2-5p and (|2-6p take the form 



"-ph,pp(a;) = -capp(a;)nph(a;) + epa(a;), ( 2-47) 

?^±,Pp(p±) = -CQ;pa(p±)n-±(p±) + epp(p±)- ( 2-48) 

By comparing with equations (|2-32p and (|2-4ip we find the absorption coefficients and emissivities to be 

app(x) = ctt Spp(x) A^ph, epa(a;) = 47r c ctt N^ N+ x^ jpa(a;), ( 2-49) 

apa(p±)=crTSpa(p±)^=F' Cpp (p± ) = 47r c (Tt iVph P± Jpp (p± ) ■ (2-50) 

2.6. Synchrotron radiation 

The kinetic equations describing synchrotron radiation need to be written in frame E, where we assume there is 
only tangled magnetic field (and no e lectric field). Using the Ein stein coefficients and the cross-sections describing 
synchrotron emission and absorption (|Ghisellini fc Sve nsson' 'l 99l[) . we ge t the collision terms for these processes in 
the electron/positron and photon equations (see also [Ochclkov et al.lll979r ): 



— [7pn±(p)] 



— [x^hphix)] 



dx / d^i 71P1 ■ 5(71 - 7 - x) {n±{pi) [1 + nph{x)] - n±{p) hpi,{x)} 

coll,syn Jo J'j ^ 

°° P P(x,j) 

dx d7i 7P ■ — ^(7-71 -2^) {"±(p) [l + ^ph(a;)] - n±(pi) nph(x)}, (2-51) 

oil X 

dj / dji jp ■ — 5(7 - 71 - 2;) {hc{p)[l + nph(a;)] - nc(pi)nph{x)} . ( 2-52) 

colLsyn Jl Jl ^ 

Here P{x,p) is the angle-integrated cyclo-synchrotron spectrum of a single electron, normalized to the electron cooling 
rate: 

P(a;,7)dx = -7s = ^^^P^ (2-53) 

O TOeC 

where Ub — B'^ /{8tt) is the magnetic energy density. One can readily verify that equations (|2-5ip conserve the total 
number of electrons and positrons, and that the total energy is conserved by equations (|2-5ip and (j2-52p . 

Under the physical conditions that we are interested in, the average energy (or momentum) of an emitted or absorbed 
photon is much lower than the energy (momentum) of the electron taking part in the process. The standard way is 
therefore to treat synchrotron processes as continuous cooling or heating for electrons and as an emission or absorption 
process for photons. 

We write the photon terms in the form 



Dhpi^{x) 



Ac es(a;) 



= -ca,(x)nph(a;) + -^^V' (2-54) 



Dt 

where as and Cg are cyclo-synchrotron absorption and emission coefficients, respectively. In the kinetic equation ()2-5p 
for the photon density nph(x) the corresponding term can be obtained by multiplying equation p-54p by SttX^^x^: 

?^ph,syn(a;) = ~cas{x)nph{x) +es{x). ( 2-55) 

The emissivity eg gives the number of photons emitted per logarithmic dimensionless energy interval dlnx, per unit 
volume and time and can be identified by comparing the corresponding terms in equations p-52p and ()2-54p : 

es{x) = T3- / Pi^^l) P^^cip) dp^ P{x,l) nc{p) d\np. ( 2-56) 
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Simil arly, by comparing the terms proportional to hph we identify the absorption coefficient (e.g. iRvbicki fc LightmanI 
[l979h : 

a.(x) = -l^ [ d^p[n,{p,) - n,{p)]P{a:,j) = -^ [ ^^P{x,^) jp dp , (2-57) 



Ancx^ 



c x^ 



dp 



where pi — 1/(7 — x)^ — 1 is the electron momentum corresponding to energy 71 = 7 — x and the second expression 
is obtained by expansion to the first order in x <C 7. 
In terms of the electron number density nc{p) the absorption coefficient takes the form: 



as{x) 



X3 



1 

Sttc x'^ 



-fP{x,-f) 



3nc{p) 



dhip 



dlnp. 



( 2-58) 



The synchrotron processes for electrons can be treated as continuous using the Fokker-Planck equation. It can be 
obtained from equation (j2-5ip employing the delta- function to take the integral over 71 and expanding 71P1 P{x,ji) 
and n±{pi) near p to the second order in the small 'parameter' x. Collecting the terms and finally integrating over 
the photon energy x we get 



§-^hpn^{p)] = -^ 



7s ipn±{p) -Hip) 7p 



dn±{p) 
97 



1 92 

2 5^ [^o(p)7P'^±(p)], 



where 



H{p) 



P{x,"f) nph(a;) x dx = -^ 
an 



P(x,7) 



nph(a;) d\nx, Ho{p) = / P{x,'^) x dx 



( 2-59) 



( 2-60) 



To get the total electron energy gain/loss rate, one has to multiply equation (|2-59p by SirX^ jdj and integrate. 
Multiplying equation (j2-54p by SttXq x^dx and integrating gives the corresponding rate for photons. Using expressions 
(l2-53p . ()2-56p . ()2-57p and ()2-60|) . we can verify that energy conservation is maintained when switching from equation 
( 2-5ip to the continuous approximation (|2-59p . 

No t e that the last t erm on the rhs of equation (|2-59p is missing in similar equations derived previously (jMcCravl 
119691 : iGhisellini et al.l [1988) . It corresponds to the diffusion due to spontaneous emission, but does not contribute 
to the electron cooling/heating. However, in most cases we expect its contribution to be negligible compared to the 
other terms. It is of the order x/j smaller than the cooling term with |7s| and, when electrons are mildly-relativistic, 
self-absorption becomes important, riph ^ 1 and the term containing H dominates. Therefore, we neglect the term 
with Hq in our simulations. Thus, for the distributions n±{p), equation p-59p takes the form 



n±,syn{p) 



d 



dlnp 



Ac,sy'a{p)n±{p) - Bc,syn{p) 



dn±{p) 



dlnp 



where 



A, 



>b) 



3^H{p) 



Bc,syn{p) = H(p)-^. 



( 2-61) 



( 2-62) 



It is worth mentioning here that other emission/ absorption processes, e.g. bremsstrahlung, can be implemented 
analogously to the synchrotron radiation, once the emissivity function of a single electron P{x, 7) (which now may 
depend on the particle distribution) is specified. 



2.7. Coulomb collisions 
The RKE accounting for electron (positron) evolution due to Coulomb scatterings is 



p-Vn±{p) 



2 [ d^pi d^p[ d^p' 






71 7i 7 
The invariant reaction rate for M0ller scattering (i.e. e 



^(^1 + £ ^ £1 ~ £') -^Coui [nc{p'x)n±{p') - nc{Pi)n±{p)] . 



and 



i^Co 



1 



'^6+) is given by (jBerestetskii et al]|1982D 
^'^^^^ +4 

(e-i)(a-i) 



( 2-63) 



( 2-64) 



and the scalar products of particles' four- momenta are defined as S, = p ■ p' and ^1 = p, ■ p' ■ As discussed bv I Baring 
(|1987f ) and ICoppi fc BlandfordI (1990), the corresponding rates for Bhabha e^e^ scattering are nearly the same in 
the small-angle scattering approximation, we therefore do not distinguish between electrons and positrons in these 
equations. 

Although the Coulomb process is coUisional in nature, it is customary to treat it in the Fokker-Planck framework, 
i.e. as a continuous diffusive energy exchange mechanism. This is due to the well-known divergent nature of the 
Coulomb cross-section for small- angle scatterings with negligible energy exchange per event, while in the parameter 
regimes we are interested in, a large number of such scatterings dominates the energy gain or loss rate of a particle 
over a much smaller number of large-angle scatterings. In frame E, where the particle distributions are approximately 
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homogeneous and isotropic, we can therefore write the Coufomb terms in the form of the Fokker-Planck equation in 
(scalar) momentum space 

d 



n±,Con\{p) 



Ac,Cou\{p)n±{p) - Bo,Coui(p)^7j 

o mp 



dlnp 
with coefficients given by 

7Coui7 d flj Dcou\ \ n / \ 1 7^^Coui 



( 2-65) 



Ac^CoAp) = ^ZJi TT 9 2 ' B^.Con\[p) = 7^ J • (2-66) 



The energy exchange rate and the diffusion coefficient can be obtained by calculating the first and second moments of 
equation (|2-63p keeping only small-angle scatterings and are expressed as integrals over the particle distributions: 



7Coui = / a(7,7i)no(pi)rflnpi, Dcoui(p) = M(7,7i) no(pi) rflnpi. (2-67) 

The rates 0(7,71) and ^(7,71) have been calculated bv lNavakshin fc Melial (|1998t l and are given in Appendix[Fl 

3. NUMERICAL TREATMENT 



We numerically solve the set of coupled integro-differential equations of the general form (|2-7[) - (|2-8p . We first define 
an equally spaced grid in the logarithms of particles' momenta: 

lnpi = lnp,„in + (j- 1) • Ap, ie[l,i,„], (3-1) 

\nxi=\nxrair. + {l-l)-^x, ^e[l,/™]. (3-2) 

Writing all differentials and integrals on the finite grids, we get three systems (for photons, electrons and positrons) 
of linear algebraic equations of the general form 

^^^^^-f:Mt:^^/-.\^^^+4), (3-3) 

i'—l 

where Ai^ is the size of the fc-th (variable) timestep. Such semi-implicit differencing scheme, where both sides of the 
equation are centered at timestep k + 1/2, is known as the Crank-Nicolson scheme (see e.g. iPress et al.|[l992) . All 
physics is contained within the matrix Mi^ii , which can be explicitly calculated at each step. The systems of equations 
for all types of particles are solved stepwise, alternating between equations and requiring a matrix inversion at every 
step. After solving a set of equations for photons, the updated photon distribution is used to calculate matrix M for 
electron and positron equations. Then we solve for distributions of electrons/positrons and substitute it to the photon 
equation and so on. 

3.1. The Chang and Cooper scheme 

The matrix Mi^ii of the linear system can be decomposed into two parts arising from the differential and integral 
terms in equations (|2-7p - (|2-8p . The differential part contributes a tridiagonal matrix, the form of the equation (e.g. 
for electrons), giving rise to it, is 



n?+^ - n? 1 



Aife Ap 

where the momentum space flux is given by 



pfe+1/2 _ j^k+1/2 

1+1/2 ^i-1/2 



3-4) 



fc+1/2 fc+1/2 



pk+l/2 _ <fe+l/2 fc+1/2 _ p,fc+l/2 '"-i+l n. 

^i+1/2 ~ 1+1/2 '1+1/2 ^i+1/2 A ■ \ ■^ ^) 

Z_ip 

The distribution function between time gridpoints is defined according to the Crank-Nicolson scheme as (omitting the 
momentum index) 

n^-+i/2 ^ i (n*^-+i + n'') . ( 3-6) 

We al so have to somehow define the distribution function between momentum gridpoints. Following I Chang &: Coopeii 
(|197Cl| ) we introduce a parameter Si so that (now omitting the time index) 

ni+i/2 = (l - Si)ni+i + SiUi, Si e [0,1]. (3-7) 

The basic idea of the Chang and Cooper scheme is to employ this parameter to ensure that the differencing scheme 
converges to the correct equilibrium solution independently of the size of the gridstep Ap . Assuming that the momen- 
tum space flux through the boundaries vanishes, the equilibrium solution tells us that it must vanish everywhere, i.e. 
F = 0. From equations P-Sp and P-7p we then have 

rij+i ^ Sj A,_n/2 Ap + -Bi+1/2 
rii Bi+i/2 - (1 - Si) Ai+i/2 Ap' 
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while the exact solution gives (jChang fc Coopeilll970D 



rii+i 

= exp 



A, 



-1/2 



1^^^ 



3-9) 



i+l/2 

We can see that using either centered-differencing {5 — 1/2) or forward differencing i5 = 0, equations (|3-8p and 
agree only to the first order in AlS.p/B. To make the correspondence exact, one has to equate the two equations and 
solve for 5i, to get 

1 1 A+l /2 

k = T^^' ^. = --^±^Ap. (3-10) 

Wi exp(wi) - 1 -Bj+i/2 

As ide from converging to th e correct equilibrium solution, such choice of 5i also guarantees positive spectra, as shown 
bv lChang fc Coopei) (|1970D . Although this method apphes to purely differential equations, we can still use it in our 
integro-differential equations to ensure that the differential part tends toward its own correct equilibrium solution, 
which would also be the correct solution for the full equation in the region where the differential terms happen to 
dominate. 

3.2. Treatment of Compton scattering 

Accurate numerical treatment of Compton scattering over a wide range of energies is not straightforward. This is 
caused by the well-known fact that at different energies the process takes place in different regimes. If the energy of 
a photon in electron rest frame is much smaller than the electron rest energy, the process takes place in the Thomson 
regime and the electron loses a very small amount of its energy in one scattering. Correspondingly, there is a sharp peak 
in the electron redistribution function R^ near p — pi. We cannot therefore numerically resolve Re on our finite grid 
and have to treat the energy loss process as continuous. On the other hand, for scattering in the Klein-Nishina regime 
the electron can lose a significant amount of its energy in one scattering. Wishing to include both regimes, we need a 
way to switch from the continuous approximation (implying a differential equation) to direct calculation of scattering 
through the integral terms. Similar treatment is required for photons, although the continuous approximation is 
only needed in the regime where the photon energy is much lower than the electron rest energy and the electron is 
non-relativistic. 

3.2.1. Scattering of electrons: separation of regimes 

Let us first look at the electron redistribution function (|2-26p . We wish to know what is the lowest incoming photon 
energy x^{pi) that can cause a shift in electron momentum pi by |Alnp|. This energy is related to the lower limit 
(I2-28|) of the integral in equation (I2-26P . If the shift is small enough, we can write 

„±/„ ^ .,^ ... ., N 1 ., iA.,i , lA„l^ ., 1„ iAi„„ I f. , Pi 



^.(Pi) « a;47,7i) = - (±|A7| + |Ap|) « -pi \A\np,\ (^1 ± ^ j > ( 3-11) 

where we have used pdp = 7^7. The plus sign applies when the electron gains energy and the minus when it loses it. 
We see that for high energy electrons, the minimum energy of photons for which we can resolve up- or downscattering 
is vastly different. However, since the upscattering (energy increase) of relativistic electrons is extremely inefficient, 
we concern ourselves only with being able to resolve their downscattering (i.e. cooling) and so use the minus sign in 
equation p-lip . Choosing jAlnpij comparable to our grid step (we use somewhat arbitrarily 3Ap) in the electron 
equation, we then state that scattering of electrons on photons with xi < x^{pi) cannot be resolved. 
We now split the redistribution function into two parts according to whether we can or cannot resolve it on our grid 

R,ip,Pi) ^R^ip,Pi) + R^ {p,Pi), ( 3-12) 

where for the first term the integral in equation (|2-26p is taken over xi < x~{pi), and the second term is defined by 
integrating over the remaining xi. To totally isolate scatterings that undergo on photons with energies below and 
above x~ , we have to write the extinction coefficient as an analogous sum, so{p) — Sq{p) + Sq{p), where 

^fT ^P)^— J ^e iPi , P) Pi dii , (3-13) 

in accordance with equation (|2-29p . For the terms containing R^ and Sq in the electron equation, we compute the 

integrals through the discrete sums, but the terms containing R^ and Sq have to be accounted for by continuous 
energy exchange terms in the equation. Since we also want to treat thermalization by Compton scattering, these 
terms have to contain a second order derivative of the electron distribution (a diffusive term). Therefore, we take the 
standard form of the Fokker-Planck equation 

N±AieMl) = -^ |7cA^±(7) - 2^ [i?e(7)A^±(7)]| , ( 3-14) 

while the exact equation for the (<) terms comes from equation (|2-25p . written here for A^±(7) 

^± coll cs(7) = -cax So<(p) A^ph^±(7) + ^^ca^N^^p f ^ X{p,Pi) A^±(7i). ( 3-15) 

J 7i 
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In order to make a physically sensible correspondence between these two representations, we demand that the first 
three moments of equations (I3-14|) and (|3-15|) were identical. Substituting equation p-13p to ()3-15|) we find 

Pi 



= 47rcf7T^phy d-fj d7i^(7l-7') R^ {pi,p) N±{^) = ca^N^^ j ^7 (tI - 7') ^^ (p) ^± (7) , (3-16) 

where similarly to the moments of the photon redistribution function (NP94), we defined the moments of the electron 
redistribution function 

^ls<ip) = — j p,-fl d7i i?e bi,p)- ( 3-17) 



The zeroth moment (giving zero in the rhs of eq. |3-16| ) is just a statement of particle number conservation, while the 
first moment gives the total rate at which the electrons gain (or lose) energy. The moments defined by equation (j3-17p 
can be calculated analytically using the exact Klein-Nishina scattering cross-section. For photons this was shown by 
NP94, while the extension of these calculations to the electrons is given in Appendix [Cl 
The moments of the continuous approximation (|3-14p are 



TV. 



±,diff,cs 



(7)^7 = 0, 



N±,diS,cs{l)ld'i= / 7cA^±(7)c?7: 



N. 



±,diff,csl. 



(7)7^^7= / [277c-f i?o(7)] N±{-f)d-f. 



( 3-18) 
( 3-19) 



( 3-20) 



Here we have assumed that the distribution function iV-|-(7) vanishes at the boundaries of integration. Exact corre- 
spondence with equation p-16p can be made if we identify 



7c = caTiVph(7i-7)sob)> i?c(7) = cfiTA^ph (71 - 7)' ^o (p), (3-21) 

while for the zeroth moment the correspondence is automatic. These moments can be computed using equations (JCIII) 
and (jC12[) . Finally, we write equation (|3-14p through n±{p) and in the form that can be included in the Chang & 
Cooper differencing scheme together with other terms 



n±,diff,cs(p) 



d 



where 



^o,cs(p) 



7c7 

9 

pZ 



d_ 



d\np 

l7^e(7) 



Ac.csip)n±(j)) - Bccsip) 



dn± (p) 



p- 



Bc.csip) 



dlnp 

2 p4 



( 3-22) 



( 3-23) 



3.2.2. Scattering of photons and three-bin approximation 

Insufficient resolution of numerical calculations can become an issue also for the scattering of photons if the electron 
energies are low enough. A photon will then exchange very little energy with an electron upon scattering and the 
redistribution function is strongly peaked near x = xi. To overcome this we propose the following approach. We 
separate scatterings that take place within some narrow interval around the energy of the incoming photon from those 
invoking photon energy outside this interval. We then approximate the scatterings taking place within the central 
interval by a continuous process and account for this by differential terms calculated through the exact moments of 
the redistribution function. 

To keep the correspondence to the electron equation, we rewrite the photon evolution equation (|2-16p in terms of 
Nphix): 



Afph,con,cs(a;) = 47r c cftN, 
dxi 



dxi 



^ ^1 

^ph(2;i) — Rph{x,xi) - A^ph(a;) — Rph{xi,x) 

Xl X 



A^ph(a;i) — i?ph(a;,a;i) - iVph(a;) — i?ph(xi,a;) 

Xl X 



( 3-24) 



where the extinction coefficient is expressed explicitly through i?pij using equation ()2-20p . Here S stands for the interval 

[xe~^'' , xe^"] and ^ means integration from to 00 excluding that interval. The width of the central region {25x in log 

units) is somewhat arbitrary, but should include at least a couple of bins, with our choice being three, i.e. 5x ~ fA^,. 

For the second integral in equation (|3-24p we wish to write a continuous approximation similar to equation (|3-14p 

d { Id 1 

iVph.diff.cs(a;) = -^ i xc Nph{x) - -— [Dph{x) Nph{x)] ^ . ( 3-25) 
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Similarly to what was done for electrons, the coefficients in equation (|3-25p are determined from the requirement 
that the first three moments of the differential and integral equations coincide. The moments of the 'central' part of 
equation p-24p (denoted by e) are 

^plcoUcsWa^'rfa; = 47rcaTiVe [ dx [ dxi {x\ - x') —Rph{xi,x)Nph{x), ( 3-26) 

' Jo Je X 

where the integration limits for x and Xi in the first term were switched, because for constant Sx the area on the 
(x, xi) plane is the same. The moments of the differential equation are similar to what were obtained for electrons 

^^ph,diff,cs(^)^^ = 0, (3-27) 

/•OO 

^ph,diff',cs(^) xdx^ X, 7Vph(a;) dx , ( 3-28) 



^ph diff cs(a^) x^dx= [2xx, + Dph{x)] Nphix) dx. ( 3-29) 



ph,diff ,cs 
"'0 

Equations p-26p and p-27p - p-29p give identical expressions for the first three moments of the 'central' part of the 
equation if we identify 

Xc = 4:TTcaTNc / dxi (xi ~ x) — R„hixi,x), Dp\Jx) = in carNc / dxi (xi — x) — Rr,h(xi,x). (3-30) 

Je X J(z X 

The 0-th moment is identically zero for both equations (|3-26p and (|3-27p . implying particle conservation. 

The integrals in equations (|3-30p are computed numerically at a finer grid. At low photon energies, the redistribution 
function can be narrower than the whole integration interval, and integration can present a problem. In this case, 
however, we can extend the integration limits in equations (J3-30I) from to cxd and to calculate the moments of the 
redistribution function analytically (NP94). Using the limits on 7^,, given by equation (|2-19p . one can show that 
scattering takes place entirely within the central interval e for incident photons and electrons satisfying the following 
relations: 

^x _ ^x 

a;<y, P<P^{x) = ^-x- (3-31) 

We can write the moments of the redistribution function in a way similar to equation (|3-17p : 

x\sq{x) = — / x^^^ dxi i?pij(xi,a;), ( 3-32) 

where the < superscript signifies that only electrons with p < p^{x) are taken into account. Equations (j3-30p then 
(for X < 5x/'2-) become 



Xc = ccttNc {xi- x) s^[x), Dp\,{x) = ccttNc {xi- xY Sq{x), (3-33) 

and can be computed using equations (|C5p - (|C6p . 
For numerical differencing equation p-25p has to be written in the form 

dnpi^{x) 



%h,diff,cs(a;) = -TTj 

omx 

where 



^ph,cs(a;)nph(a;) - i?ph,cs(a;)- 



9 In a; 



( 3-34) 



^-<^'^|-SG^)' W..)^i^^ (3-35) 

3.3. Pair production and annihilation 

The numerical treatment of pair-production and annihilation processes in our code is fairly straightforward. The 
only potential difhculty can arise from the non-linearity of the absorption term in the photon equation. To deal with 
this we have chosen the simplest possible approach: for calculating the pair-production opacity at each step we simply 
use the photon distribution from the previous step. The error caused by doing so is not expected to be significant in 
most cases. It is well-known that a photon of energy x will most efficiently interact with photons of energy xi « 3/a;, 
thus if its energy is not very close to the electron/positron rest energy, the photon will most likely annihilate on another 
photon of a vastly different energy than its own. Therefore, we can visualize two separate populations of photons that 
pair-produce on each other, with the dividing energy at rUcC^. The photon distribution from the previous step is then 
taken to be the 'target' population on which the photons that are being evolved pair produce. 

Since we wish the numerical scheme to treat electrons and positrons identically (particularly when we are dealing 
with pure pair plasma), while at each step one of them has to be evolved first when the outcome of the other is yet 
unknown, we use a fully implicit scheme for the pair annihilation terms. 

The quantities Rj~^{j^,x,Xi) crpa(7+,7-) and crpp(a;,a;i) defined by equations (|2-37p . (|2-39p and (|2-46p are precal- 
culated on a fine grid and thereafter averaged within the electron/positron and photon bins used by the code. The 
integrals in the expressions (|2-36p . (|2-38p . (|2-44p and (|2-45p for emissivities and absorption coefficients are calculated 
through discrete sums. 
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3.4. Treatment of synchrotron processes 

One of the main difhcuhies in numerically treating synchrotron processes in compact sources is that the optical 
thickness of the medium due to self-absorption might become extremely large at low energies compared to, say, 
Thomson optical thickness. Almost all photons that are produced are immediately absorbed, so very few escape. But 
the energy which those few carry away comes from the small net energy exchange rate between electrons and photons, 
which we need to keep track of to maintain the energy balance. Near the equilibrium, in the photon equation we have 
two large terms describing emission and absorption, which nearly exactly cancel out. A small error in either of them 
produces a significant error in the total energy transfer rate. In the electron equation this transfer rate is given by the 
difference in the synchrotron cooling and heating rates. To maintain the energy balance between the two equations, 
we need to ensure that in our numerical scheme these rates are seen identically by both equations. 

In discretized form, the synchrotron processes for electrons/positrons are described by equations (|3-4p - (|3-5p . with 
n = n±, A = Ae,syn and B = _Be,syn- To obtain the total energy gain we have to multiply equation (j3-4[) by 7i Ap, 
sum over i and sum the corresponding terms in the electron and positron equations. Assuming vanishing boundary 
currents, we have 
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^i+1/2 "o,i+l/2 ~ J^i+1/2- 
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( 3-36) 



where A7i+i/2 = 7i+i — 7i and we have omitted the time index k + 1/2 for brevity. The exchange rate as seen by the 
photon equation can be evaluated by writing the integrals in emissivity and absorptivity expressions (j2-56[) and (|2-58p 
as sums over the grid, multiplying equation p-55p by xi Ax and summing over I: 
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where Pi^i — P{xi,pi). Changing the order of summation, identifying the sum over the photon distribution as the 
discretized version of the definition H{p), and noticing that J2i Pi.iXiAx gives the electron cooling rate — 7s,i, we get: 



A^, 



ph 



AtkAV 



= Ea. 



7s 



37»g» 



liHi »c. i+i - ric, i 
p1 Ap 



( 3-38) 



To make equations p-36p and p-38p identical (except for the sign) we have to make subtle changes in the definition 
of coefhcients and the way integrals are numerically calculated. In equation p-38p we have to define the coefficients in 
between the electron momentum gridpoints, at i -t- 1/2, substitute the electron distribution rio ^ by n^^ i+1/2 (except in 
the derivative term), where the latter is calculated using the same Chang & Cooper coefficients Si as in the electron 
equation, and sum up to z = i^ — 1 instead of z„i. This amounts to defining the emission and absorption coefficients 
as 
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( 3-39) 



Also, the coefhcients A and B entering the momentum space flux 
rate P-36P should be written as 



5]) and thus also the electron energy exchange 
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( 3-40) 



1/2 A7,+i/2 vr / i+1/2 

and ensure that the energy exchange rates as seen by the electron 



which become identical to (|2-62p in the limit Ap 
and photon equations are the same. 
The only discrepancy left is that we cannot use the same n , and ric in both equations. This is because each 

of them contains a function n'^+^, which, in the equation that we evolve before, is not known for the other type of 
particle. The solution to this, at least in the average sense, is to regard the time-grids for each equation as shifted by 
a half timestep. Then n^^^ obtained from one equation can be used as n^+^i'^ in the other and vice versa. 



3.5. Coulomb collisions 

Coulomb scattering only redistributes the energy between different parts of the lepton population. It is easy to see 
that the total energy is conserved in the sum of two equations (|2-65p for electrons and positrons, provided that 0(7, 71) 
is antisymmetric, the latter simply reflects the energy conservation in two-body interactions. Similarly to synchrotron, 
our numerical treatment has to ensure that the conservation is exact, otherwise unphysical runaways can occur near 
the equilibrium. 
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The flux in momentum space in equation (|3-4p for Coulomb scattering is given by equation (j3-5p with coefficients 
expressed as (see eq. |2-66p 
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( 3-41) 
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The total energy exchange rate is identical to equation (|3-36p for synchrotron. 
Let us now look separately at terms containing 7 and D. For 7 we have 
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It is now easy to see that this quantity can be made to vanish if we write 
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provided that a is antisymmetric. The terms containing -D(7) in the energy exchange rate are 
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where we have redefined the coefficient B as 
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( 3-43) 



( 3-44) 



( 3-45) 
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One can see that equation (|3-44p has the form of an integral over a full differential and, as such, should vanish 
provided that D = at the boundaries. To ensure this numerically for any electron distribution we write explicitly 
'^6,2+1/2 = (1 ~ Si)nc^i+i + SiUe^i and demand that the coefficient in front of ric^i in equation (j3-44p is equal to zero for 
every i. Rearranging terms, we get 
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The expression in the curly brackets is identically zero if we set 



P / i+1/2 
while the boundary terms S~ and S^ vanish if 
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( 3-46) 



( 3-47) 



( 3-48) 



Using expressions p-43p in the first term in coefficient A and equations ()3-47p and (j3-48p in the definition ()3-45p . we 
ensure precise energy conservation in the numerical scheme. 

4. NUMERICAL RESULTS 

Our careful treatment of the micro-physical processes makes the code applicable over a wide range of parameter 
regimes. The current version covers 15 orders of magnitude in photon energy (from 10"^ to 10^" eV) and 8 orders 
of magnitude in electron momentum, while there is no fundamental difficulty in extending this range further, e.g. to 
TeV energies for application to blazars. Energy conservation is achieved to within 1% in the majority of cases. All 
the rates and cross-sections of different processes have been precalculated once and for all and are read into memory 
as the code initializes. A typical simulation for 200 grid points in photon energy and electron momentum on a 3 GHz 
PC running Linux takes between a few minutes and half an hour. 

In order to test the performance of our code in different parameter regimes, we have chosen three setups from earlier 
works and run the code with similar parameters for comparison. 
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Fig. 1. — Equilibrium (a) photon spectra and (b) electron distributions (Thomson optical depth per Inp, i.e. ne(p)a"T^) for various 
stochastic heating compactnesses Jtij as labeled. The size of the emission region is ij = 10^"* cm, the soft input radiation has a compactness 
Is = 10 and a blackbody temperature Tbb = 15 eV, the injectio n compactness is /nth = 10- The thin solid line on the right panel shows a 
Maxwellian fit of temperature Te = 53 keV. Compare to fig. 1 in lCoDpH ||1999I V 

4.1. Non-thermal pair model 

As a first test we compare our code to the well-known pair plasma code eqpair bv lCoppil (| 19921 [19991 ). eqpair also 
considers an uniform emission region into which high-energy electrons/pairs are injected, mimicking an unspecified 
acceleration mechanism. Some low-energy photons are also injected, emulating a source of external soft radiation (e.g. 
accretion disk). The high-energy pairs cool by Compton scattering and Coulomb energy exchange with colder thermal 
pairs. The Compton upscattered photons can produce electron-positron pairs which then upscatter more photons etc., 
initiating a pair cascade. Once the pairs cool down to low enough energies, the timescale of the systematic energy losses 
becomes longer than that of diffusive processes, leading to relaxation into a low-energy thermal distribution. In eqpair, 
Coulomb collisions between particles are assumed to be the thermalizing mechanism. However, the thermalization 
process is not treated entirely consistently in this code in a sense that there exists only one thermal bin into which 
particles are put once they have cooled below a certain threshold energy, chosen to be 7 = 1.3. The electron temperature 
associated with this thermal bin is nevertheless calculated self-consistently from energetic considerations. Furthermore, 
the code does not consi der thermalization by synchrotron self-absorption, which can be an efficient mechanism if the 
medium is magnetized (jGhisellini et al.lll988l GHS98). 

The setup of this test run is similar to what was used for fig. 1 in ICoppil (|1999( ). We switched off synchrotron 
processes in our code and left other processes. We inject a Gaussian distribution of pairs centered at 7inj = 10'^ and a 
low-energy blackbody distribution of photons. In addition, there is a background electron plasma present with optical 
depth Tp = 0.1. There is no escape term for pairs, meaning that all injected pairs eventually annihilate transferring 
their energy to the radiation field. The power injected as non-thermal pairs is parametrized by compactness 



(Tx L 



nth 



'nth 
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(4-1) 



where Lnth is the injected luminosity (including rest mass) and R is the linear dimension of the emission region. 
Similarly, we define the compactness of the injected soft radiation as 

is 
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(7X 



.c^ R 



(4-2) 



where L^ is the relevant luminosity. To mimic acceleration with less than 100% effi ciency, a dditio nal power is supplied 
to low-energy electrons in the form of continuous heating, parametrized by /th- In iCiJPpi! (|1999f ) this energy was just 
given to the thermal bin, but since we do not have such bin in our code, we need to explicitly specify the form of this 
heating. This is done by stochastic acceleration prescription of the form 
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(4-3) 



The momentum diffusion coefEci ent is assumed to t ake the form characteristic of stochastic acceleration by resonant 
interactions with plasma waves (jDermer et al.lll996[) . £'acc(p) oc p'^. We have chosen q = 2 in our calculations. The 
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Fig. 2. — Evolving (a) photon spectra and (fe) electron distributions (t(p) = ct'^ Rna(p) / p) for Gaussian electron injection under action of 
Compton and synchrotron processes at different times (in R/c units) as labeled. The source size is i? = 10^^ cm, the magnetic compactness 
is ^B = 10 Enid the injection compactness i„tii = 1. Compare to fig. 1 in GHS98. 



mean energy gain rate of a particle resulting from equation (|4-3p is 

dt I 



stoch. 



^^ Wd^^M] , 



4-4) 



where j3 = p/7 is the particle speed. We can see that for a power-law diffusion coefficient the gain rate is proportional 
to p"^"^ in the relativistic regime, while in the nonrelativistic regime it is proportional to p'. Our choice q — 2 means 
that at high energies Compton losses always overcome gains by stochastic acceleration, the main effect of the latter 
process is therefore the heating of low-energy pairs. 

The differential term given by equation (j4-3p is included in the Chang & Cooper scheme on the same grounds with 
other continuous terms. Therefore before discretization it has to be written in the form compatible with equations 
and (13^: 
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(4-5) 



The results of the test are shown in Fig.[TJ Varying the amount of stochastic heating ( /th) kee ping all other parameters 
co nstant, we see that we can well reproduce the behavior of the spectrum in fig. 1 in iCoppj |l999). Just as expected 
bv iCoppil (|1999f ). the equilibrium electron distribution is hybrid: Maxwellian at low energies with a nonthermal high- 
energy tail. Note that we get such shape even if we switch off Coulomb scattering. The thermal-looking distribution 
is produced by the stochastic heating itself, which gives a Maxwellian slope at low energies irrespective of the shape of 
Dacc{p), while the location of the peak of the distribution is determined by the balance between heating and Compton 
cooling. The behavi or of the spe ctrum in response to varying the power of stochastic heating seen in Fig. [TJa was 
analyzed in detail bv lCoppil (|1999D . we are not going to repeat it here. 

4.2. Thermalization by synchrotron self-absorption 

For the second test, we compared our results with these of GHS98. They studied electron thermalization by syn- 
chrotron self-absorption in the presence of Compton cooling. The electron cooling, heating and diffusion due to the 
synchrotron were described by equation (|2-59p (without the last term), while Compton scattering was assumed to take 
place in the Thomson regime and contribute only to systematic cooling. Furthermore, the treatment was not fully 
self-consistent since only the electron equation was actually solved. While the equilibrium synchrotron spectrum was 
self-consistently calculated at each timestep from the formal solution of the radiative transfer equation, the Comp- 
tonized spectrum was not. Thus only the synchrotron spectrum entered the electron heating rate by self-absorption, 
while the radiation energy density needed to account for Compton cooling was estimated from energetic considerations. 

We ran our code with the same parameters used to obtain the results in figs. 1 and 2 in GHS98. The pair 
production/annihilation and Coulomb scattering have been switched off for this test. High-energy electrons are injected 
into the emission region, with the total power (including rest mass) parametrized by the injection compactness ^nth- 
The magnetic compactness is defined by 

Ib = ^^RUb, ( 4-6) 



Radiative processes in hot magnetized plasma 



17 



' L ' ' ' I ' ' ' I ' 
'.. \ 100 




-6 -4 -2 

log {x=E/m^c ) 



Fig. 3. — Equilibrium (a) plioton spectra and (6) electron distributions for injection ||4-7I| for various injection compactnesses i^th as 
labeled. Parameters: R = 10^^ cm, (b = 30. Compare to fig. 2 in GHS98. 

where Ub is the magnetic energy density. In addition there is an external source of soft blackbody photons assumed 
to arise from reprocessing half of the hard radiation by cold matter in the vicinity of the emission region. The electron 
escape timescale is fixed at tcsc = R/ c. 

In the first case the injected electrons have a Gaussian distribution peaking at 7 = 10. The evolution of this 
distribution is followed in time as it cools and thermalizes by Compton and synchrotron processes. We can see that 
our results shown in Fig. [2] are almost identical to those presented in fig. 1 in GHS98. However, we would like to stress 
that we also compute self-consistently the photon spectrum. We see the partially self-absorbed synchrotron bump at 
small energies, then the blackbody photons and two Compton scattering orders at higher energies. 

In the second case we calculated the steady-state particle distributions for different injection compactnesses. The 
injected electron distribution (per unit Inp) is 



Qc = Qo ^ exp 1 , 



7 



7c 



(4-7) 



where 7c — 3.33. The resulting equilibrium electron distributions plotted in Fig. [3^ are again very similar to the ones 
obtained by GHS98 in their fig. 2. The corresponding radiation spectra shown in Fig.[31a are computed self-consistently 
and simultaneously with the electron distribution (while the spectra in fig. 4 of GHS98 are calculated a posteriori, 
i.e. after the equilibrium electron distribution has been determined). As discussed in GHS98, if the source is strongly 
magnetically dominated, the equilibrium distribution is almost purely Maxwellian. When the injection compactness 
increases, Compton losses become non-negligible and the electrons cool down to lower energies before they have time 
to thermalize. Notice that at the highest compactness (^nth = 100) the temperature of the Maxwellian part of the 
distribution inferred from Fig. [3^ deviates appreciably from the one obtained by GHS98. This is caused by the fact 
that at high compactness a significant fraction of the soft radiation is Compton upscattered to energies comparable to 
the energies of the Maxwellian electrons. These photons are therefore not effective in cooling the electrons any further. 
However, in GHS98 Compton cooling is accounted for through a term proportional to the radiation energy density, 
which includes all photons, and therefore overestimates the cooling rate. Overall, the simple prescription for Compton 
cooling without actually solving the photon equation appears to work well in the parameter regimes considered here. 

4.3. Gamma-ray bursts from stochastically heated pairs 

Finally, we compare our code to the Large Particle Mo nte Carlo c ode by 'Ster n et al.l (|1995( ) , with all the processes 
operating now. The setup is similar to the one used in St ern fc Pou tancn (2004) for simulating the spectral evolution 
of gamma-ray bursts. They consider an initially optically thin distribution of electrons in a cylinder-shaped emission 
region. Arguing that impulsive first-order Fermi acceleration would result in cooling spectra that are too soft to be 
consistent with observations, energy is instead supplied to the electrons continuously, mimicking dissipation by plasma 
instabilities behind the shock front. As electrons are heated to relativistic energies in the prescribed background 
magnetic field, they emit synchrotron radiation, providing seed photons for Compton upscattering. The high-energy 
upscattered photons then initiate pair-production. 

In our simulation we consider a spherical region permeated by magnetic field and start by heating a cold electron 
distribution (with initial Thomson optical depth tq = 6 x 10^^) according to the stochastic acceleration prescription 
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Fig. 4. — Evolving (a) photon spectra and (fe) Thomson optical depth per Inp for stochastically heated pairs at different times (in 
units R/c) as labeled. Parameters: the source size R = 10^^ cm, the magnetic compactness Ib = 0.3, the stochastic heating compactness 
Ith = 30, the initial Thomson optical depth of electrons is tq = 6 X 10~^. For t = 0.1,0.3 we also plot positrons, at later times only the 
electrons as their opacities are nearly identical. Compare to fig. 2 in lStern &: PoutanenI II2004I ). 

(j4-5|) . No pair escape is allowed. The results of simulations are shown in Fig. 3] and can be compared to a similar 
iig. 2 in lStern fc Poutanen (2004). In both cases the electrons are rapidly heated to about 7 ^ 100, as determined by 
the balance between stochastic heating and synchrotron cooling. As the photon field builds up, additional cooling by 
Compton scattering causes the electron 'temperature' to start dropping. After about 1/3 of the hght crossing time, 
the number of photons upscattered to the MeV range becomes large enough to start significant pair-production. With 
the increasing pair density (at i = 1, opacity has grown by a factor of 20) the available energy per particle decreases, 
causing a further drop in the temperature of the now almost pure pair plasma. After about ten light-crossing times 
the Thomson opacity is tt = 1.3 and the pair density reaches the value where the pair annihilation and creation rates 
are balanced and a steady state is attained. 

The spectral behavior seen in Fig. [4la is similar to what was obtained by lStern fc PoutanenI (|2004f ). The synchrotron 
peak rises first, being initially in the optically thin regime and thus following the evolution of the peak of the electron 
distribution according to a; oc 7^. The first Compton scattering order lags slightly behind synchrotron, while the second 
scattering order is initially in Klein-Nishina regime and thus hardly visible at all. As the electron temperature drops 
and the peak of the first scattering order evolves to lower energies, the second order shifts to the Thomson regime and 
becomes comparable to and eventually dominant over the first order. At the same time the decreasing temperature 
and increasing pair opacity causes the synchrotron emission to switch to optically thick regime and the synchrotron 
luminosity to drop dramatically. The plasma becomes photon starved and the Comptonized spectrum hardens. 

5. CONCLUSIONS 

We have developed a new computer code for simulations of the radiative processes in magnetized rarefied plasmas 
encountered in the vicinities of accreting black holes and relativistic jets in active galaxies and gamma-ray bursts. We 
take into account Compton scattering, pair production and annihilation, synchrotron processes and Coulomb scattering 
without any limitations on the energies of the photons and electrons/positrons. We solve coupled integro-differential 
kinetic equations describing time evolution of the photon and electron/positron distributions simultaneously. The 
equations contain both integral and up to second order differential terms. The Fokker-Planck differential terms are 
substituted when necessary instead of the integral terms with coefficients computed exactly from the moments of 
the integral equation. This allows us to study thermalization of the lepton distribution by Compton and Coulomb 
scattering and synchrotron self-absorption. Processes involving bremsstrahlung can be easily added to the code, while 
for the conditions considered in the paper they are not important. 

The presented technique guarantees energy (and particle, when relevant) conservation with high accuracy which 
is especially important when dealing with strongly self-absorbed synchrotron radiation. The implementation of the 
Chang and Cooper scheme gives the correct shape of the particle distribution at low energies. The area of application 
of the code is enormous as it can deal with photons and leptons covering many orders of magnitude in momentum 
space, with no potential difficulty of extending it to even lower/higher energies. We present a number of test runs, 
where we consider problems previously solved by other methods. We compute non-thermal pair cascades, and study 
lepton thermalization by synchrotron self-absorption, as well as model the emission from the stochastically heated pairs 
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that might have a relation to the prompt emission of gamma-ray bursts. We find a good agreement in the parameter 
space where comparison is feasible while the differences can be explained by our improved treatment of microphysics. 
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APPENDIX 

A. RELATION BETWEEN THE COMPTON REDISTRIBUTION FUNCTIONS FOR PHOTONS AND 

ELECTRONS 



The redistribution functions defined by equations (|2-18p and (|2-27p can be written as 



i?ph(a;,a;i,7i) = -7-^ Pi / p d-y d^fl 6{'^i + xi - 7 - a;) / d^fii d^wi F S{pi + Xi - p - x) 
i?e(7,7i,a^i) = -7-^ xi j X dx d^uj (5(7i + xi — 7 — x) j d^ili d^uji F 5{p^ + Xi — p — x) 



(Al) 



(A2) 



We see that the inner integrals are identical in both expressions. Because of rotational symmetry, the only angle left in 
the calculation after performing the integrals over d^flid^uji is the angle between the momenta of outgoing particles. 
Therefore we can write d^fl = d'^u = 27Td(, where C = fl-Lj. We also see that d7(5(7i-|-xi— 7— x) = dx6{'yi+xi—j—x), 
so we find from equations (lAip and (|A2p that the redistribution functions are related as 

ppiRe{j,ji,xi) =xxiRphix,xi,ji), (A3) 

where one of the energies/momenta has to be replaced from the condition a; + 7 = xi + 71. 

B. COMPTON REDISTRIBUTION FUNCTION 

The isotropic Compton redistribution function defined in equation (j2-18p can be written as an integral over the 
scattering angle (NP94) 

_ /•/i+ '^+ 

i?ph(a;,a;i,7i) = / R{x,Xi,ji,fi)dfi = T{x,Xi,ji,^j.) 



The limits of integration are given by 

{—1 if |a; — xil > 2xx\^ 
— 1 if jx — xil < 2xx\ and 71 > 7*(x, xi, —1), 
/^_ if \x- xx\ < 2xxi and 7m < 7i < l*{x,xi,-l), 

where 
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pi +li{xi - x) +PiV(7i +xi - xY - 1, 7m = 1 + {x- Xi + \x-xi\)l2. 



The quantity ^^^{x^xi, —1) is the minimum electron energy needed to scatter a photon backwards (i.e. /i 
xi to x: 



(Bl) 
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(B3) 
-1) from 

(B4) 



7*(x,a;i, -1) = [x — xi + {x + xi)^/l + llxxi\/2. 

The angle- depend ent r edistribution function j ?(a;, x i, 71, ii) was first derived bv lAharonian fc AtovanI (|198lD. se e also 
iPrasad et al.l (|1986l ) and lNagirner fc PoutanenI (|1993l ). The angle-averaged function was obtained bv lJoneJ ([Till), but 
the presented expressions are very cumbersome a nd the loss of accu racy occurs for small photon energies and large 
electron energies. An alternative function given by iBrinkmannI (|1984f) and NP94 does not suffer from these problems. 
We use here the latter expressions. The primitive function T(a;, a;i,7i, /x) can be expressed through functions of one 
argument as 
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where w = 1 — fi and Q ~ \/{x — xi)^ + 2xxiw. The functions H are given by the differences 

H = A{h^) - A{h+), Hn = An{h^) - Anih+), 

where 

A{h)^VTTh, h+^[i-fi+xi)^-l]w/2, h^ = [{ji-xy~l]w/2. 

The zeroth function Aq is 

°^ ' \arcsin(x/^)/V^ if /i < 0, 



(B5) 
(B6) 

(B7) 

(B8) 
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while the others can be expressed through its derivatives as 



A„{h)^{-2r^^^4-\h), (B9) 



and can be computed by the recurrent relation 

1 r 9n -I- 1 1 1 

(BIO) 



An+l{h) = - 



2n + l ^ ,^, 1 

\2n-l\ "^' A^"+^ih) 



or for \h\ < 1 via series 

A (h) 1^"'^' V (2n + 2fc-l)!! {-h^ 

^"^''^ - (2n-l)!! ^ ^^)n 2n + 2fc + l- ^^^^^ 



fc=0 



Direct computations using (jB6|) lead to numerical cancellations at small photon energies x,xi ^ 1. NP94 describe in 
details how they should be dealt with. 

C. MOMENTS OF THE COMPTON REDISTRIBUTION FUNCTION 

The moments of the photon redistribution function given by equation (j3-32p can be written explicitly using equa- 
tion (P^H]) as 

—-, ^ 3 21 fd^pd^piffxi ., ~.xp.4 fn^\ 

a;lso a;) = TTrTTl^- / x[nc{p)FS^, CI 

16tt A^iVc a; 7 7 71 xi 

where we have denoted S^ — 6{p + x^ — p — x) for brevity. We now define (NP94) 

(A) soiO = ^U^^ A FS^ (C2) 

lOTT £, J Xl 71 

and 

M^, 7) = ^^^^ J d'n e (xl) so(0, (C3) 



where ^ — p ■ x. Using equations (jC2p and (IC3I1 . we get for equation (ICl 

2 



a;lso(x) =47r^^3-^x* / il^{p) pUp ^ ,{x , -f) . (C4) 

Analytical expressions for ^^ along with asymptotic formulae for different limiting cases can be found in NP94. For 
calculating x^ and -Dph(x) using equations p-33p we need terms like {xi — x)* so(x), for « = 1,2, which are simply 

(xi-a;)so(a;)=47r^3-^ x / no(p)p2 dp (*i - *o), (C5) 



(xi-x)2so(a;)=47r-^^— a;2 / ne(p)p' dp (*2 - 2*i + *o). (C6) 

We now rewrite the moments of the electron redistribution function given by equation p-17p using equation ()2-23p 

—- f . 3 2 I [ (fix (Pxi (Ppi , ^ / N p e4 tn-7\ 

167rA£7Vph7J x Xl 71 
We can define quantities analogous to equations (|C2[) and (jC3[) : 



3 1 f (fxi cPpi 



(7l)so(0 = t|-j j'^^^^^\F5- (C8) 

167r^J Xl -^ 



(C9) 



71 
and 

*^(^>7) = ^^^/rf^-e(7l)5o(0- 
Equation (|C7P then takes the form 

7Jso(p) = 47r ^3^ 7' / nph(x)x2dx $,;(x,7), (CIO) 

while the terms needed for calculating -^^ and D^ (7) using equation (|3-2ip are 

(7i-7)so(p)=47r ^3^ 7/ nph(x)x2dx($i-$o), (Cll) 

2 



(7i-7)2so(p)=47r^^3-^7^ / nph(x) x^ dx ($2 - 2$i + $o)- (C12) 



Radiative processes in hot magnetized plasma 

From considerations of energy conservation one would expect a relation between the rat es 
(IC12|) . To see this, consider the quantities x (^i — ^o) and 7 ($1 — $0) entering equations (jC5|) and (IClip 

1 
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x{^i 



*n 



d-n ^(xi -a;)so(0, 



where 



Analogously for electrons 



Aii^x 

Sir (Pxi cPpi 



3 1 [d^x 
{xi -x)sQ[i) = — -- / 

IbTT^J Xl 



71 



(a;i ^x)F5'^ 



7($i-$o) 



1 



d^uj ^(71 -7)so(0> 



47ra;7 

3 1 f d^xi d^pi 



, ,,, 3 1 f d-'x 

(7i-7>^o(e) = Y^^y — 



71 



(71 - 7) F 54. 



C5|). ea and (pTj) . 
(C13) 

(C14) 

(C15) 
(C16) 



Due to the energy conservation (5- function, a;i — a; = 7 — 71, and we thus get 

(a;i-x)so(e) = -(7i-7)so(C)- (C17) 

Also, after performing the integrals over d^xi d^pi in equations (|C14p and (jC16[) the only remaining angle that 
{xi — x)so{^) and (71 — 7)so(^) can depend on is the one between the incoming photon and electron momenta. 
Therefore in equations (|C13[) and (|C15[) we can write d'^fl — d'^uj — 2TTdC,, C, ^ Vl ■ u and conclude that 



(C18) 



x(*i-*o) = -7(*i -*o)- 
Using the same arguments one can show that 

x^ (*2 - 2*1 + *o) - 7^ (*2 - 2$i + $o). (C19) 

We can thus use the analytic expressions for '^i for calculating the rates 7c and I?c(7) for electrons as well as photons. 
Note that since '^q = $0 we also have analytic expressions for calculating the total scattering cross-section for electrons 
through equation (|C10|) . setting i — Q. 

D. ELECTRON-POSITRON PAIR-PRODUCTION AND PAIR-ANNIHILATION RATES 

The quantity R^~f entering both the pair-production and the annihilation rate expressions (J2-36I) and p-44p can be 
written as 

1 r n TT^ TT- 1 ™ 

+ T(7_,a:,a:i,Xcm) + r(7_,xi,a;,a;cm) ■ (Dl) 



i?^T-(7_,X,Xl) 



'^J{x + xiY 






The primitive functions are (jNagirner &: Loskutovl[T999f) 

T(7_,a;,a;i,a;cm) 



(a;a;i)3/2 h 






2(a;a;i)3/2 



X {xi + x) + 7- {xi — x) ~ 2x1 



'ixxiAa{h) 



(D2) 



where h — [(7- — x)"^ — 1] x^^/xxi and a;cm is the photons' energy in the center of momentum frame. Upon exchanging 
X and Xl in the preceding expression one also has to reverse the arguments in function h. The quantities A and ^o 
are identical to the functions defined by equations (|B7p an d (IB8D for Compton scattering. The expressions similar to 
(jDip and (|D2p have been derived bv iSvenssonI ()1982[ ) and iBoettcher fc Schlickeiserl (|1997|) . However, their formulae 
suffer from cancellations when h approaches zero, while in equation (|D2p cancellation appears only in the term 
[ylo(/i) — A{h)]/h, which can easily b e co mputed via Taylor series for small h. 
The integration limits in equation (jDip are 



■^cm -^cm and X^jjj 



minl-^xxi, x^^^ 



where 



(■^cmj ~ (,7cmj r, 



7-7+ + 1 ± v/(7--l)(7i-l) 



(D3) 
(D4) 



and 7+ = a: 4- Xl — 7- (which follows from the energy conservation). 
The lower limit of integration in equation (|2-36p are expressed as 



/{[2x-7-(l+/9-)]7-(l + /3-)} ifx>x+. 



j(^) = i7_(l-/3_), xf' = L/{[2a;-7-(l-/3-)]7-(l-/3-)} if ^ < 



(D5) 



l-x 



if x_ < X < x+, 
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where we have defined 



x± = -[l + 7-(l±/3-)] 



(D6) 



Because we always have x\ > x^^', the latter sets a lower limit for the energy of either photon for producing an 

electron (or positron) of energy 7_. 

The lower limits of the momentum integrals in equation (I2-44|) are given as follows () S venssoni 1 1 9821 ) : 

'7^-) if a; < 1/2, 
(L) _ J 7^"-' if 1/2 < x < 1 and 7+ < 7b, (l) _ /7A if 2; < 1/2, 



7^'*') if a; > 1 and 7+ < 7b, 
1 in all other cases. 



where 



v(±) - 



F± 



1 

F^ 



F±=2a;-7+(l±/3+), 7a 



7+ = 



4x" 



if a;> 1/2, 



4x 



7b 



2x2 _ 2a; + 1 
2a; -1 



The total pair-production cross-section (in units of ctt) is given by (jGould fc Schredeiill967HZdziarskilll988[ ) 



tTpp(a;,a;i) 
where 



3 1 






In w ^ — ^, In^ w + 2 In (w + 1) + 4Li2 



v + 1 



y/v + l 
xxi — 1 and w 



vr 
"3" 



^/v + l + s/v 



and Li2 is the dilogarithm defined by 



Li2(r) = - f 
Jo 



V^u + T 

ln(l - s 



ds. 



The total pair annihilation cross-section is fomid to be (e.g. ISvenssoiilll982r ) 



crpa(7+,7-) = o 



1 



/?L7cm^(/3cm) - 272^ + -L^p,^) 



(D7) 

(D8) 

(D9) 
(DIO) 
(Dll) 

(D12) 



8 7+7_p+p_ 
Here we have defined L{P) — ln[(l -I- /3)/(l — /3)], while the limits of integration are given by equation (|D4p and 

E. CYCLO-SYNCHROTRON EMISSIVITIES 

The cyclo-synchrotron emissivity (here in units s^^ str"-'^) at photon energy x in t he directi on given by angle 9 to 
the magnetic field for an electron moving at a pitch-angle a with velocity /3 = p/7 is (jPacholcz vk 1970) 



'q{x,0,a) = -— ttf x^ V 



cos 9 ~ /3 cos a 
sin6' 



Jf(z)+/32sin2a J/^(z) 



5 N a;[l — /3 cos a cos( 

7 



(El) 



where ai — e^ /ch is the fine-structure constant, b = B/Bcr is magnetic field in units of the critical field Ba = 
mlc^ /{eh) = 4.41 x 10^^ G, J; and J[ are the Bessel function and its derivative, and their argument z = xp sin a sin 9 /h. 
Averaging over pitch-angle and integrating over 9, we get the angle-averaged cyclo-synchrotron spectrum 



P{xn)^7i 



dcosa 27r 



dcos9 ri{x, 9, a). 



(E2) 



Direct summation over harmonics works fine for mildly relativistic electrons 7 < 3. In this c ase, we first use the S- 
f unct ion to integrate over the energy bin, and then integrate numerically over the angles (see e.g. lMarcowith fc Malzad 
l2003f l and sum over harmonics contributing to a given bin. The same procedure is used for any larger 7 at photon 
energies x corresponding to the first 30 harmonics (i.e. a; < 30 6/7). At higher x, we use two different methods. In the 
ultra. - relativistic regime 7 > 10 we use the angle-averaged relativistic synchrotron spectrum (Crusius fc Schlickeisen 
[l98fiHclfisellini et al.iri988h : 



P(x,7) 



SV^cttUb 1 



rrirC b 



T X < K^/3{x)Ki/3{x) - -X K^/^{x) - K^/^{x) 



(E3) 



where x = a;/ (37^6) and Ky is the modified Bessel function. For 3 < 7 < 10, we substitute the sum over harmonic in 
equation (|Eip by the integral over I and use the ^-function t o take it. The angu l ar inte grals are then taken numerically 
Alternatively we use the approximate formulae proposed bv lKatarzvhski et al.l ([2006') . which ignore harmonics. These 
give identical results for the simulations presented in the paper, because low harmonics are self-absorbed. 
All the emissivities P(x,"f) are renormalized to guarantee the correct cooling rate given by equation (|2-53[) . 
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F. COULOMB EXCHANGE RATES 
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The Fokker-P l anck treatment of Coulomb (M0ller) scattering in relativistic plasma was first considered by 
iDermer fc Liana (|1989f) . Useful analytical expressions for the energy exchange rate and diffusion coefficient for (small 
an gle) scattering of tes t elect ron of energy 7 interacting with the background electrons of energy 71 have been derived 
bv lNavakshin fc Melial (|1998f ). These expressions (corrected for a few misprints and reorganized) are: 



and 



SccTTlnA 

a(7,7i) = T 71 -7 X 7,71 

4 771PP1 

,, , SccTTlnA 

«(7,7i) = 7 ^(7,71) 



4 771PP1 



where In A '--^ 20 is the Coulomb logarithm, 



X(7,7i) 



7rcl 



Trol 



Prol (7rcl - 1) 



d7, 



rcl 



Prcl 



7rcl + 1 
Prcl 



+ ln(7rel + Prcl) 



Trol 



and 



/•7+01 72 

A(7, 71) = / -^ (7rtl - 7rcl)(7rcl - 7rel) '^7rcl 



.-, P 



(Fl) 
(F2) 

(F3) 

(F4) 



rcl 



= ^ - ( 7' + 7? + ^ j ln(7rcl + Prcl) + —^ [7rcl(7' + ll) " 2771] + p,,i (2771 - ^ 



Here Pioi = ■\/7rci ~ ^ i^ ^^e relative momentum and the integration limits are 

7rei=77i±PPi- (F5) 

We can obtain simple approximate expressions for the energy exchange and diffusion coefficients by approximating 
the integrals in equation (jFSp using one-point trapezoidal and in equation (JF4| a 3-point Simpson's rule: 



f N_3 , (71 " 7) 771 

•2(7: 71) ~ nCO'T In A = 

2 (771 - l)v(77i)^ - 1 



and 



d(7, 7i) « C(Tt In a ■ 



ihp^pI 



[{ill? - if^ ' 
which agree with the exact coefficients reasonably well, except in the region 7 « 71. 



(F6) 
(F7) 
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